# Extreme weather events do not increase political parties' environmental attention


#'   This script
#'    > performs PanelMatch models based on dictionary classification
#'    > creates one graph for all models
#'    > exports Extended Data Figure 2


# Tim Wappenhans, António Valentim, Heike Klüver, and Lukas F. Stoetzer
# First: 22-02-2023
# Last: 10-11-2023


# PACKAGES --------------------------------------------------------------------
if (!require("pacman")) install.packages("pacman")
pacman::p_load(
  tidyverse,
  PanelMatch,
  panelView,
  ggpubr,
  ggrepel,
  cowplot,
  here
)


# DATA WRANGLING FOR PANELMATCH ------------------------------------------------
df.panelmatch <- read_csv(here("data/data_output/press_weekly.csv"))  |> 
  mutate(
    
    # only integers allowed 
    parlgov_id = as.integer(parlgov_id),
    t = as.integer(t),
    country = unclass(as.factor(country_name)),
    family = unclass(as.factor(family_name_short))
    
  ) |> 
  
  select(
    
    # Covs
    parlgov_id, t, family, country, cabinet_party,
    
    # Treatments
    treat, treat_flood, treat_storm, treat_fire, treat_temp,
    
    # Outcome
    climate_prs, pr_total
    
  ) |> 
  
  mutate(
    
    # climate
    outcome_climate = climate_prs / pr_total,
    outcome_climate = ifelse(is.nan(outcome_climate), 0, outcome_climate),
    
  ) |> 
  
  as.data.frame() |> 
  
  filter(!is.na(parlgov_id)) 


# Globally define leads and lags, iterations
lag <- 6
lead <- 6
n_iterations <- 1000


# create empty tibble for data viz
gg.did <- data.frame(
  coef = numeric(), 
  se = numeric(),
  model = numeric()
)


# LOOPS TO CREATE ONE GRAPH ACROSS ALL MODELS ----------------------------------
models <- c("Overall",
            "Green",
            "Left",
            "Right",
            "Opposition",
            "Government",
            "Floods",
            "Storms",
            "Wildfires",
            "Extreme Temp.")

for (model in models) {
  if (model == "Overall") {
    df.tmp <- df.panelmatch
  } else if (model == "Green") {
    df.tmp <- df.panelmatch |> filter(family == 5) 
  } else if (model == "Left") {
    df.tmp <- df.panelmatch |> filter(family %in% c(3,6,8)) # com, lib, soc
  } else if (model == "Right") {
    df.tmp <- df.panelmatch |> filter(family %in% c(1,2,4,7)) # agr, christian, Cons, Radical Right
  } else if (model == "Opposition") {
    df.tmp <- df.panelmatch |> filter(cabinet_party == 0)
  } else if (model == "Government") {
    df.tmp <- df.panelmatch |> filter(cabinet_party == 1)
  } else if (model == "Floods") {
  df.tmp <- df.panelmatch |> mutate(treat = treat_flood)
  } else if (model == "Storms") {
    df.tmp <- df.panelmatch |> mutate(treat = treat_storm)
  } else if (model == "Wildfires") {
    df.tmp <- df.panelmatch |> mutate(treat = treat_fire)
  } else if (model == "Extreme Temp.") {
    df.tmp <- df.panelmatch |> mutate(treat = treat_temp)
  } 
  
  PM.results <- PanelMatch(lag = lag, 
                           lead = 0:lead, 
                           data = df.tmp,
                           time.id = "t", 
                           unit.id = "parlgov_id",
                           treatment = "treat", 
                           outcome.var = "outcome_climate",
                           match.missing = TRUE,
                           refinement.method = "mahalanobis",
                           covs.formula = ~ I(lag(outcome_climate, 1:lead)),
                           exact.match.variables = "family",
                           size.match = 10, 
                           qoi = "att",
                           forbid.treatment.reversal = FALSE)
  
  
  ## Treatment
  PE.treat <- PanelEstimate(sets = PM.results, 
                            data = df.tmp,
                            se.method = "bootstrap",
                            number.iterations = n_iterations)
  
  # Store Coefs and SEs
  gg.treat <- summary(PE.treat, 
                      verbose = F, 
                      bias.corrected = F) |> 
    as.data.frame() |> 
    mutate(model = model,
           period = row_number() - 1) |> 
    rename(low_ci = `2.5%`,
           up_ci = `97.5%`) |> 
    select(-std.error)
  
  ## Placebo
  PE.placebo <- placebo_test(PM.results, 
                             data = df.tmp, 
                             number.iterations = n_iterations)
  
  # placebo estimates  
  placebo_estimates <- PE.placebo$estimates |> as_tibble() |> 
    mutate(period = row_number() - (lead + 1)) |> 
    rename(estimate = value) 
  
  
  # placebo CIs, percentiles from bootstrapping
  percentiles <- PE.placebo$bootstrapped.estimates |> as_tibble() |> 
    summarize(across(everything(), ~ quantile(., probs = c(0.025, 0.975)))) 
  
  percentiles_longer <- percentiles |> 
    mutate(quantile = row_number(),
           quantile = case_when(quantile == 1 ~ "low_ci",
                                quantile == 2 ~ "up_ci")) |> 
    pivot_longer(cols = starts_with("t-"), names_to = "period", values_to = "value") %>%
    pivot_wider(names_from = "quantile", values_from = "value") |> 
    mutate(period = substr(period,2,3),
           period = as.numeric(period))
  
  # join placebo estimates and placebo CIs
  gg.placebo <- left_join(placebo_estimates, percentiles_longer,
                          by = "period") |> 
    mutate(model = model)
  
  # bind placebo to treatment
  gg.tmp <- rbind(gg.placebo, gg.treat)
  
  # Add to main viz DF
  gg.did <- rbind(gg.did, gg.tmp)
  
}


# PLOTS ------------------------------------------------------------------------
# globally define settings
pre_color <- "#4885C1"
post_color <- "#AE3A4E"

pre_fill <- "#4885C1"
post_fill <- "#AE3A4E"


pre_shape <- 17
post_shape <- 16

# create indicator
gg.did <- gg.did |> 
  mutate(pre = ifelse(period < 0, "pre", "post"))

# reorder
gg.did$model <- factor(gg.did$model, levels = models) 



# Store it so don't have to bootstrap every time anew
write_csv(gg.did,
          file = here("data/data_output/ggdid_dic.csv"))

# pick back up
gg.did <- read_csv(here("data/data_output/ggdid_dic.csv"))

# data wrangling for plot
gg.did <- gg.did |> 
  mutate(
    
    # make it percentage points
    estimate = estimate * 100, 
    low_ci = low_ci * 100, 
    up_ci = up_ci * 100, 
    
    # add label for ggplot
    lab_text = case_when(
      model == "Overall" & period == 0 ~ 
        "No increase in attention\n towards climate after\n extreme weather events"
      )
    )

# Store different model types for graphs
## Models
overall <- c("Overall")
families <- c("Green", "Left", "Right")
govt_oppo <- c("Opposition","Government")
disaster_types <- c("Floods", "Storms", "Wildfires", "Extreme Temp.")

figmodels <- list(overall, families, govt_oppo, disaster_types)

## names for plots
gg.names <- c("overall", "families", "govt_oppo", "disaster_types")

## titles within plots
gg.titles <- c("", "Party Families", "Government Control", "Event Type")

# Loop all plots
i <- 1
for (figmodel in figmodels){
  gg.did |> 
    filter(model %in% figmodel) |> 
    ggplot() +
    
    # Gray pre treatment shade
    geom_rect(aes(xmin = -Inf, xmax = -0.5, ymin = -Inf, ymax = Inf),
              fill = "#e0e0e0", alpha = 0.2) +
    
    # CIs
    geom_linerange(aes(x = period, ymin = low_ci, ymax = up_ci, color = pre),
                   alpha = .7,
                   linewidth = 0.75) + 
    
    # Coef
    geom_point(aes(x = period, y = estimate, shape = pre, color = pre),
               size = 3) + 
    
    # fixed lines
    geom_hline(yintercept = 0, linetype = 'dotted') + 
    
    # look
    ylab("Percentage change in attention") +
    xlab("Weeks relative to extreme weather event") + 
    scale_x_continuous(breaks = gg.did$period, labels = gg.did$period) +
    facet_wrap(~ model, 
               scales = "free", 
               ncol = 2) +
    scale_color_manual(values = c("pre" = pre_color,
                                  "post" = post_color)) +
    scale_shape_manual(values = c("pre" = pre_shape, 
                                  "post" = post_shape)) +
    scale_fill_manual(values = c("pre" = pre_fill, 
                                 "post" = post_fill)) +
    guides(shape = "none",
           fill = "none",
           color = "none") +
    
    coord_cartesian(xlim = c(min(gg.did$period), max(gg.did$period))) + 
    
    # theme
    theme_bw() +
    theme(panel.grid = element_blank(),
          plot.title = element_text(hjust = 0.5)) +
    ggtitle(gg.titles[i]) ->  gg.tmp
  
  assign(paste0("fig.", gg.names[i]), gg.tmp)
  
  i <- i + 1
}



# Annotate Overall graph individually
gg.did |> 
  filter(model %in% overall) |> 
  ggplot() +
  
  # Gray pre treatment shade
  geom_rect(aes(xmin = -Inf, xmax = -0.5, ymin = -Inf, ymax = Inf),
            fill = "#e0e0e0", alpha = 0.2) +
  
  # CIs
  geom_linerange(aes(x = period, ymin = low_ci, ymax = up_ci, color = pre),
                 linewidth = 1.15) + 
  
  # Coef
  geom_point(aes(x = period, y = estimate, shape = pre, color = pre),
             size = 5) + 
  
  # fixed lines
  geom_hline(yintercept = 0, linetype = 'dotted') + 
  
  # look
  ylab(" ") +
  xlab(" ") + 
  scale_x_continuous(breaks = gg.did$period, labels = gg.did$period) +
  facet_wrap(~ model) +
  scale_color_manual(values = c("pre" = pre_color,
                                "post" = post_color)) +
  scale_shape_manual(values = c("pre" = pre_shape, 
                                "post" = post_shape)) +
  scale_fill_manual(values = c("pre" = pre_fill, 
                               "post" = post_fill)) +
  guides(shape = "none",
         fill = "none",
         color = "none") +
  
  coord_cartesian(xlim = c(min(gg.did$period), max(gg.did$period))) + 
  
  # theme
  theme_bw() +
  theme(panel.grid = element_blank(),
        plot.title = element_text(hjust = 0.5)) ->  fig.overall


# Annotation for first graph
fig.overall_an <- fig.overall + 
  #  Annotation
  geom_curve(
    data = gg.did |> filter(model == "Overall", period == 0),
    aes(
      x = 0.95,
      xend = period + 0.15, # nudge
      y = 2.25, # nudge,
      yend = 0 + 0.5 # nudge
    ),
    color = "gray",
    alpha = 1,
    curvature = 0.1,
    size = 0.75,
    arrow = arrow(
      length = unit(0.025, "npc")
    )
  ) +
  # for each annotation, add a label
  geom_text(
    data = gg.did |> filter(model == "Overall", period == 0),
    aes(
      x = 1.5,
      y = 3.5,
      label = lab_text,
      color = "gray",
      hjust = 0.5,
      vjust = 1
    )) 
  
fig.overall_an

# bring all together, keep only one x-title
cowplot::plot_grid(fig.overall_an + theme(axis.title.y = element_blank(),
                                          axis.title.x = element_blank()),
                   fig.families + theme(axis.title.y = element_blank(),
                                        axis.title.x = element_blank()),
                   fig.govt_oppo + theme(axis.title.y = element_blank(),
                                         axis.title.x = element_blank()),
                   fig.disaster_types + theme(axis.title.y = element_blank() ),
                   ncol = 1,
                   rel_heights = c(2, 2, 1.1, 2))  +
  # manually add y-title
  draw_label("Change in issue attention to climate", 
             size = 12, angle = 90,
             x = 0) +
  # increase margins so it's completely visible
  theme(
    plot.margin = margin(t = 0, r = 1, b = 0, l = 1, "cm")) -> fig.all


# EXPORT -----------------------------------------------------------------------
fig.all 
ggsave(here("results/graphs/fig_ex_data_02.pdf"),
       width = 7,
       height = 15)
